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1.  INTRODUCTION 


In  this  report  we  give  a numerical  procedure  for  integrating  the  bivariate  normal  density  func- 
tion over  a convex  polygon.  The  Fortran  IV  computer  program*  developed  from  it  is  fast  and  is 
designed  to  yield  the  output  probability  to  3,  6,  or  9 decimal  digit  accuracy.  As  far  as  we  know,  it  is 
the  fastest  and  most  versatile  program  of  its  kind  - most  versatile  in  the  sense  that  it  handles,  with 
three  prespecified  levels  of  accuracy  in  the  output,  arbitrary  convex  polygons**  rather  than  just  tri- 
angles and  quadrilaterals.  We  make  note  at  this  time  that  the  program  serves  as  a basic  subroutine  for 
the  automatic  computation  of  the  bivariate  normal  over  an  arbitrary  polygon.  A complete  program 
for  this  much  more  general  case  has  been  written,  checked  out,  and  is  operational.  Us  description  is 
deferred  to  a later  report. 

Our  procedure  for  the  convex  polygon  case  depends  on  a fast  method,  with  prespecified  accu- 
racy, to  evaluate  the  bivariate  normal  distribution  over  an  angular  region,  A.  In  particular,  we  wish 
to  evaluate 


U) 


P(A)  = 


(1  ~P2)~ 
2 v Oj02 


. (w -/iJU-Oj) 
2P -z-z 


where  (p, . p2)  is  the  mean  and 


lpo,o3 


PU,Oj' 

(t 2 

J 


the  covariance  matrix  of  the  normal  random 


variable  (\v,2)  with  correlation  coefficient  p.  The  angular  region  A is  defined  as  the  semi-infinite 
part  of  the  plane  bounded  by  two  intersecting  directed  straight  tines.  Of  course  by  this  definition 
there  are  four  such  regions,  and  therefore  it  is  necessary  to  always  state  which  of  them  is  involved. 


The  well-known  linear  transformation 


(2) 


X a 


Ft5 


y a 


* - H 


reduces  the  integrand  of  ( 1 ) to  one  with  circular  symmetry,  namely 
(3)  I* (A)  “ P(A)  *-v4Scxp|~<xJ  + yJ>/2J  dx  dy, 

***  A 

*THe  program  i*  coded  for  the  CDC-6700.  a large-scale  binary  computer  capable  of  one  million  operations  per 
second,  li  has  a 60  bit  binary  word  length  of  which  48  arc  used  to  express  the  mantissa  of  a number, 

•M  ite  term  convex  polygon  will  always  mew  a closed  convex  polygon. 
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where  A,  like  A,  is  an  angular  region,  since  (2)  takes  straight  lines  into  straight  lines.  Thus  we  deal 
only  with  (3)  hereafter  unless  noted  otherwise. 

An  extensive  literature  exists  on  methods  for  integrating  the  bivariate  normal  variate  over 
various  simple  geometries,  where  the  ultimate  objective  is  to  appropriately  utilize  these  integrations 
to  evaluate  the  distribution  over  a polygon,  [ 2,3 ,4, 5,6,7 ] . One  such  case  is  where  A in  (1)  forms  a 
right  angle  at  (h,k)  with  the  sides  of  the  angle  directed  parallel  to  the  w and  z directions.  When  the 
mean  is  zero  and  the  variances  are  equal  to  one,  ( 1 ) for  the  angular  region  just  described  is  denoted 
by  <l>(h,k,p)  and  is  called  the  bivariate  normal  integral  [3]  or  the  bivariate  normal  probability  func- 
tion [9,  p.  936] . We  shall  make  reference  to  «t>  in  Section  6.  We  show  it  is  equivnlent  to  (3)  where 
the  given  right  angle  is  transformed  to  an  angular  region  A and  then  show  that  a recent  method  for 
computing  «J>,  (3] , is  slower  than  our  procedure  for  obtaining  the  same  result  from  (3). 

The  idea  of  integrating  over  an  angular  region  seems  to  have  originated  with  Gideon  and  Gur- 
land  (hereafter  G & G),  (4]  * [5] . As  observed  by  them,  the  idea  of  integrating  over  an  angular 
region,  as  expressed  by  (3),  is  a natural  and  easily  visualized  way  to  obtain  the  probability  over  a 
polygon.  In  Section  5 we  shall  discuss  and  compare  their  computing  method  with  ours. 

In  Section  3 we  show  how,  by  utilizing  (3)  over  a set  of  angular  regions,  we  obtain  the  proba- 
bility over  a convex  polygon.  Our  approach  differs  here  also  from  what  G & G advocate.  In  Section 
7,  we  give  some  numerical  results.  The  computer  program  is  described  in  Section  4 with  its  Fortran 
listing  given  in  Appendix  1).  In  the  next  section  we  give  some  analysis  and  also  the  algorithm  for 
evaluating  (3).  Its  implementation  into  a computer  program  is  not  straightforward  since  certain 
precautions  are  necessary  as  will  be  explained  in  Section  4. 


2.  ALGORITHM  FOR  P(A(R,  0 , , 6>3 )] 

In  this  section  we  derive  the  algorithm  by  which  we  evaluate  (3);  i.e.,  we  obtain  that  part  of  the 
circular  normal  distribution  over  the  angular  region  A(R,  0, , 0,)  as  the  shaded  region  shown  in 
Figure  1.  Lines  (T)  and  (D  form  the  boundaries  of  this  region.  R denotes  the  distance  from  the 
origin  to  the  vertex  of  A(R,  0, , 02). 

It  is  convenient  because  of  circular  symmetry  in  the  integrand  of  (3),  to  perform  a rotation  of  axes 
such  that  the  line  L and  the  x axis  coincide  with  A rigidly  rotated  as  shown  in  Figure  2.  Hereafter 
we  shall  always  assume  such  a rotation,  through  the  angle  C\  has  been  carried  out. 

The  coordinate  transformation 

(4)  x * R ♦ rcosd,  y a rsinfl,  101  < », 


is  used  in  (3)  to  obtain 


•Wc  are  grateful  to  Pete  Shugatt  at  White  Sands  Missile  Range,  New  Mexico  for  bringing  their  Wisconsin  report  (4) , 
to  our  attention. 


Figure  I.  Angular  Region,  A(R,  Gj , d2),  (shaded  region) 
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where 

(6)  p = R cos  0. 

An  integration  by  parts  on  the  integral  in  r yields 


(7) 
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where 


<8)  erfc(x)  ® l - erf(x) 


t)  a f ziOdt, 


J(x)  ^ o 


Using  (7)  in  ($)*  carrying  out  the  obvious  part  of  the  0 integration,  (5)  becomes 


(9) 


HA)  62  2 


*RJ/2 


2* 


f 


u|erfs»!u)/4u>)  *10  f , 


where 


(10)  u » pA/2  « (R V2)eos0. 

We  note  for  R * 0,  (V)  gives  the  exact  result  directly, 

(U)  P{A(O,A0)1  » A0/2*.  MMj  -0,. 

liquation  (9)  gives  the  relation  for  HA)  upon  which  our  program  is  based.  A similar  relation 
was originally  derived  by  Amos,  (I  J,  in  an  entirely  different  way. 

Tire  difficulty  in  evaluating  the  Integral  in  (9)  is  resolved  by  obtaining,  for  a given  6 > 0,  the 
minimax  polynomial  fit  to  etfc  (u)/?.(u)  fo.  0 < u < C(6).  Namely,  a set  of  constants  ^ and  a 
kart  positive  integer  K are  found  such  that 


4 


(12) 


erfc(u)  - z(u)  Z atuk  <-7: 6,  0 <u<C(6). 

0 * 1 \/» 

The  constant  C is  chosen,  once  5 is  specified,  such  that 

(13)  2^  $5exp[~j(x2  + y2)]dxdy  = “■  erfc (\LlyfT)  - e = 6V v, 

A 

with 

(14)  C = Rly/I,  A = A(R,-~,i), 

For  6 £ 5(~4),  5( -7),  5(-lG),  2(- 13)  the  ak  and  C are  given  in  Appendix  A.  For  example 
for  6 as  5< -7),  (=5  X I0'7),  we  give  C = 3.5505,  K = 9.  The  way  in  which  e was  chosen  in  (13) 
is  explained  below. 


front  ( 1 2)  and  (9)  we  have 


(15) 
where 

(16) 


f0,  - 0,  K / "I 

KA)  - __ \~Y—  - £ \ (yrj  I «***'«  >»j . 


■ 


j a r*4n  P*  eosk0d0, 


so  that 
(17) 

with 

(IS) 

Hence 

(19) 


J0  “ 0j  - 0,  , J,  a sin  &2  --j^-sind,. 


where  it  is  emphasized  that  f IS)  and  (19)  hold  only  when  10.1  < 9/2,  i = 1.  2.  This  folios  s from 
( 1 0),  because  u > 0 in  ( 1 2)  and  R > 0 in  ( 1 0)  imply  cos  0 > 0.  For  cases  outside  this  range,  wc 
make  use  of  the  fact  that 


5 


(20) 


P[A(R,O,0)]  = erfc  (-^  sin  0)  - P[A(R,  O,tt-0)],  ~ < 0 < * , 


where  we  prefer  to  work  with  the  coerror  function,  erfc  (see  (8)),  instead  of  the  univariate  cumula- 
tive distribution  function  of  a normal  variable.  They  are  related  by 


y erfc  (xA/2)  = 1 - J dt  • 


The  implementation  of  (19)  and  (20)  is  discussed  in  Section  4,  which  deals  with  the  computer 
program. 

2 

We  now  show  that  if  a maximum  error  of  —pr6  is  made  in  approximating 

(22)  f(u)  = erfc  (u) , 0 < u < C (8)  f 

as  noted  in  ( 1 2),  then  the  truncation  error  in  computing  P(A),  using  (19),  can  be  no  larger  than  8/^/F. 
Indeed,  from  (12) 

(23)  f F(u)  - 2 akuk  I < 5 eu\  u > 0,  F(u)  ■ M , 

o z(u) 

and  since  u > 0,  we  have 


(24)  -~e-R2/2  f9ueul  d0  < '2  f [uF(u)  - I akuk+1 1 d0  < -e'R2/2  f°2 ue“2  d0  . 

11  Jo.  * Je.  o K » Je 


But,  with  (10), 


and  (24)  then  implies 


Rj/2j(BjueU2  d0  = J^9a(^cos0)exPl~f‘sin2  <W 

= [erf(-^  sin62)  - erf  (^=  sind,))  < s/ir, 

c-r2/2  r o k 

~ Jo  luF(u)  - 2 akuk+l ) d0  < hjy/v  (=  e). 


This  accounts  for  the  way  e was  chosen  in  (13), 

The  dominant  part  of  the  computation  in  evaluating  P(A)  from  (19)  is  the  generation  of  the 
sum  of  terms  { ak  Jk  + , | . Two  situations  can  occur  for  which  this  sum  docs  not  contribute  to  the 
value  of  P(A)  to  w**hin  the  accuracy  specified,  namely  when  R is  ‘'small”  and  when  R is  ‘‘large." 
In  the  first  case,  wt  have 


•'-oxu  ».  tkiMMMMinna  iv  i>Mhi4M> * » >>t-a  j*.* 


I 

f. 


(27) 


P(A)  as  (1  - R2/2) 


I*  - 1 r 

L 2?r  ix  JQ 


where  with  uF(u)  = g(u), 
R 


, R2  2 ciy/n 

+ -x-  COS'4  6 


g(u)  = cose  (l  + -y- cos2  6)^ 

=^cos#[^^cos9  + 0<R5)] 


JgdO  de 


R 


+ 0(R  - j 


1 - ->  ~cosd  + 0(R3) 


Carrying  out  the  0 integration  in  (27),  we  obtain 


(28) 


P(A)  a f “ 2v7j  sin ^ ^ ” e>) * - 4“ 2p.)  + °«J>- 


Thus,  when 

(29) 

then 

(30) 


1 

2^ 


R . „ 

^2smd2 


-—sis-  Q 


R 1 
"/2vA 


e,  (e  s 6/y/?,  see  (13)), 


P(A)  is  A0/2v. 


Extending  the  above  analysis  one  an  show  that  the  R3  tenn  in  (28)  is  given  by 

J 


1 

v* 


sin2  0 cos  0 ( 


and  upon  integration  yields 

(31)  E H +V?  (^)  (sta>  “ Sta’®l> 

Hence  P(A)  is  approximated  to  within  e by  (28),  without  the  0(R3)  tenn,  when 

(32) 


3,/ff 


In  the  other  circumstance,  when  R is  sufficiently  large,  a parameter!?  an  be  determined, 
depending  on  e,  sudi  that  if 


(33) 


R > R (or  R2/2  > R2/2) . 


■ ‘ • V fit--2  V S* 

■ , , •-in.rt  I.  i,  ■ I « t1 
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then  P[  A(R,  0t , 02  ] < e for  I 0j  I , I d2  I < 7r/2.  So  in  this  case  that  part  of  the  computation  for 
P(A)  which  uses  (19)  can  be  omitted,  but  one  erfc  function  is  still  required  for  each  1 6. 1 > nil 
(i=  1,2)  (see  (20)).  1 


We  note  from  (13)  and  the  fact  that  P[A(R,  — 7r/2,  0)]  is  an  increasing  function  of  0,  that  for 
R > R 

(34)  P(A(R, 0j , 02)j  < p[a'R,~ I|)J  =lerfc(R/v/2),  I0j I,  !«,!< ir/2,  0,<0r 
Consequently,  we  choose  R such  that 

(35)  y erfc  (R/\/2)  ~ e - 6 fy/W  , 


and  observe  that  C from  ( 14)  and  K/y/2  are  the  same  for  a given  e.  Geometrically  it  means  that  the 
region  to  the  right  of  the  vertical  line  x = R /^/Tdoes  not  contribute  to  P(A)  to  within  the  specified 
accuracy. 


3.  USE  OF  P(A)  TO  COMPUTE  P(H) 

In  this  section  we  show  how.  using  probabilities  over  angular  regions,  the  probability,  P(H), 
over  a convex  polygon  H is  obtained.  In  (41  they  propose  using  probabilities  over  triangles  and 
quadrilaterals  to  obtain  the  same.  Our  procedure,  however,  is,  in  general  more  efficient,  as  shown 
below,  we  require  only  N angular  regions  for  an  N-sided  convex  polygon,  whereas  they  need  at 
least  3(N  - 2)  regions  if  H is  decomposed  into  triangles.  If,  for  N even,  H is  decomposed  into  quadri- 
laterals, or  quadrilaterals  plus  one  triangle  for  N odd,  then  one  needs  2N  - 4 or  2N  - 3 angular  regions, 
respectively.  We  remind  the  reader  that  our  ultimate  purpose  in  developing  a program  for  computing 
P(H)  is  tc  use  It  as  a subroutine  to  evaluate  the  probability  over  an  arbitrary  polygon.  As  stated 
earlier  this  has  been  done  and  will  be  discussed  together  with  a computing  program  in  a later  report. 

Let  M(N,  t, , . . . , t^)  denote  a convex  polygon  of  N sides  with  vertices  at  coordinate  points 

t, tN , where  tk  = (xk , yf; ) and  the  points  j t.  J are  given  in  counterclockwise  order;  i.c,  so 

that  the  area  of  K is  on  the  left  as  one  traverses  the  boundary  continuously,*  Then 

(36)  P(H)  * I*( A, ) -NS  PtAJ  + IHA*), 

where  using  Figure  3 with  N = 6,  A,  is  the  angular  region  determined  by  any  interior  angle  of 
H with  its  vertex  assigned  as  t, , A,,  i - 2, . . . , N - 1 , are  angular  regions  determined  by  the  exterior 
angles  of  H at  vertices  t2 , ... , tN . , , respectively,  as  shown  in  Figure  3 and  AN  is  the  angular  region 
obtained  from  the  vertical  angle  of  the  interior  angle  of  H at  tN . It  is  easy  to  argue  the  validity  of 
(36)  by  noting,  c.g  in  Figure  3 (N  - 6),  that  the  probabilities  over  the  disjoint  shaded  regions  E,, 
i “ 2, 3, . . . , N - 1 , excessively  diminish  the  result  for  P(H)  by  an  amount  exactly  compensated 
for  by  the  addition  of  P(AN ).  A formal  proof  of  (36)  is  not  given  in  this  report. 


•Note  that  (2)  imps  convex  polygons  into  convex  polygons. 
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Figure  3.  Convex  Polygon,  N ® 6,  Showing  Angular  Regions 


In  keeping  with  our  efforts  to  maintain  an  efficient  program,  as  described  in  the  previous 
section,  we  make  use  of  the  following  result.  Let  0(i)  denote  the  quantity  -0V  which  appears  in 
(9),  for  the  flh  regular  region  probability,  ?(Aj).  Then,  since  the  interior  angles  of  an  N -convex 
po»ygon  add  up  to  (N-2)s  radians,  have 

(V)  Om  a -0(1)  + N£  0(i). 

t-2 

Thus  our  program  generally  will  require  only  N-l  cails  to  the  tan'1  routine  instead  of  N by  using 
(37).  The  accumulation  of  the  rigid-hand  side  is  denoted  in  the  flow  charts  as  o.g.,  see  boxes 
(2,8,1 1,24,33]  in  flow  chart  U,  page  1 1 


4.  COMPUTER  PROGRAM  FOR  P(H)  (AND  P(A)) 

In  this  sect' on  we  discuss  the  Fortran  iV  program  for  the  evaluation  of  P(H)  or  P( A),  Uie 
normal  probability  distribution  over  a convex  ptdygon  H of  N sides,  or  an  angular  region  A.  res|>ec- 
tivcly.  Two  flow  charts  I,  U are  given  on  the  nex<  ’wo  pages  which  show  ihe  flow  and  major  steps 
of  the  program.  It  will  be  helpful  to  refer  to  these  charts  during  the  discussion.  A location  in  the 
flow  charts  will  be  identified  by  chart  number  and  b~x  number,  eg.  (I,iOj  refers  to  chart  1 l>  x 10. 

One  input  to  the  program  is  the  sequenc.{tk}  , where  tk  denotes  the  (x,y)  coordinates  of 
the  kth  vertex  of  H with  the  vertices  ordered  in  a counterclockwise  direction.  The  value  of  N is 
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FLOW  CHART  i 


•Use  4-cjuadrant  tan'1  touiinc  (-a,  r 


specified  with  N set  to  one  if  P(A)  is  desired.  In  this  case  3 points  are  required,  as  in  the  case  of  a 
triangle  where  N = 3,  in  counterclockwise  order,  i.e.,  so  that  the  region  A is  to  the  left  as  one 
tranverses  the  boundary  lines  with  the  only  vertex  at  tj . A parameter  is  set  specifying  whether  3,  6, 
or  9 decimal  digits  are  desired  in  the  output  P(H)  or  P(A).*  The  associated  values  of  various  param- 
eters are  given  in  Appendix  A,  namely  a;  , a[ , a2,  a3,  R/^/T.  Also  listed  in  that  Appendix  are  values 
of  these  parameters  for  P(H)  or  P(A)  computable  to  twelve  decimal  digits.  These  however  are  not 
incorporated  into  the  program  but  could  be  with  no  difficulty  if  desired. 


It  is  imperative  for  the  program  to  operate  properly  that  the  tk  be  given  in  counterclockwise 
order;  i.e.,  with  the  area  on  the  left  as  one  travels  along  the  boundary  of  H or  A.  Two  typical 
examples  are  shown  in  Figure  4,  where  P is  wanted  over  the  shaded  or  hatched  regions. 


N = 6 N =>  1 


Figure  4.  Typical  Regions  for  H and  A 

Point  t,  for  H can  be  taken  initially  as  any  vertex  point,  however  when  this  program  is  used  for 
arbitrary  polygons,  it  will  cycle  the  points  and  renumber  them  so  that  the  new  t , is  the  point  with 
the  property 

(38)  t,  = y^ , < yk  with  xJ  < xk  if  ^ » yk , k a 1, . . . , Nj  , 

(This  feature  is  not  shown  on  the  flow  charts,). 

For  A,  t,  must  be  specified  as  the  vortex  point,  as  shown  in  Figure  4 above. 

In  order  to  evaluate  P(H),  N angular  regions  |Akj  are  treated,  one  at  each  vertex  of  H,  and 
their  probabilities  |l>(Ak  )j  combined  appropriately  as  explained  in  Section  3 (see  (36)).  For  a 
particular  Ak  - A(R,dj,02),  the  inequality  below 


•We  make  note  of  the  fact  here  that  the  specified  number  of  correct  decimal  digits  In  computing  P(H)  may  not  be 
achieved  in  the  unlikely  case  that  the  errors  associated  with  a majority  of  the  angular  regions  have  the  same 
sign  and  thus  add  to  a total  error  of  as  much  as  NS. 
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(39) 


B = R2/2  < at  = ire2 


is  tested  where  is  taken  from  (29).  If  it  holds  then  [I,  5]  is  used  to  evaluate  P(A),  which  is 
then  stored  in  I, 


P(A)  = (02  - 0j  )/2ir  = ^-tair'C- 


[uz  - vwl 
uw  + vz  2 ’ 


where  the  tan"1  is  obtained  from  a four  quadrant  subroutine  which  gives  its  output  in  (-ir,ir] . The 
quantities  u,  v,  w,  z are  initially  defined  in  {I,  1 J and  subsequently  in  [II,  20,  25,  37]  depending  on 
which  angular  region  is  involved.  The  angles  02  and  6\  are  measured  in  radians  and  are  as  shown  in 
Figures  l or  2,  page  3,  with  A0  always  positive  from  0,  counterclockwise  to  02. 


If  the  inequality  in  (39)  is  not  true,  then  a rotation  of  axes  is  carried  out,  [I,  9] , as  indicated 
in  Figures  1 and  2.  Quantities  gj  /Dj , hj  /D{ , g2/D2,  h2/D2  are  computed,  [I,  10] , where 


(40) 


with 


gj/Dj  = cos 0 j -*■  g,,  h(/Dj  sindy  -*  h1  , 

g2/D2  =-~cos02  -*■  g2,  h2/D2  =~sin02  -*  h2. 


(41) 


D,  a [*2(w2  + z2)]*  , D2  = [2(u2  + v2)]*  . 


We  have  for  the  first  of  (40) 
(42)  R 


cos 


0 (4  * v 

1 ' “ Mx.  vv  + y.  z)2  + (x.  z - 


,2^ 

(xkw  + ykz)2  + (xkz  - ykw)‘] 
° (xkw  + ykz)/[2(w2  + z2)]*  « g,/Dj . 


2,W 


The  other  relations  in  (40)  arc  found  in  the  same  way.  The  location  denoted  in  the  charts  as  P 
contains  the  output  I*(A)  if  N = 1 , or  P(H)  if  N > 3.  Location  <p  contains  0 if  N s 1.  If  N > 3 
and  (33)  docs  not  hold,  then  4>  contains  0(N),  (See  (37)),  at  exit. 


In  [I,  28] , the  inequality 


(43)  B < o2  s (9*c2)1/,t  (B  “ R2/2), 

is  tested.  If  it  holds  then  P(Ak ) is  given  by  (28),  [I,  5, 32J , where  o2  is  taken  from  (32). 

In  general,  the  program  distinguishes  12  different  types  of  angular  regions  which  arc  ex- 
haustive and  arc  characterized  by  the  signs  of  the  numbers  g, , g. , h, , l>2  as  computed  in  [ 1 , 10) . 
Examples  of  each  of  the  1 2 regions  arc  shown  below  in  Figure  5 with  the  terms  used  to  evaluate  P, 


• **"..*:  . -v  ' - :*v  ■ aresveaKSi 
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(P  = P(A)).  It  is  assumed  the  rotation,  [I,  9]  has  been  carried  out  as  described  above,  so  that  the 
vertex  of  A is  on  the  positive  x-axis  (not  at  the  origin).  The  angle  between  the  directed  lines  labelled 
(T)and  (2)  is  always  measured  from  (T)  to  (2)  in  the  counterclockwise  direction  and  it  is  non-negative 
and  always  no  larger  than  tt  (sin  (02  - 0 1 ) > 0),  since  we  are  dealing  with  convex  polygons.  We  allow 
it  <2n  only  if  N = 1 . In  this  case  we  evaluate  P(E-A)  = P,  where  E denotes  the  entire  plane, 
and  find  P(A)  = 1 - P.  The  boxes  that  apply  for  N = 1 only,  showing  the  details  just  mentioned,  are 
[I,  36,  37,  38] . In  the  situation  shown  in  Figure  5,  we  denote  the  probability  over  the  angular 
region  between  (T)  and  (2)  by  P,  and  note  in  the  expressions  for  P below  each  diagram,  that  if  g;  < 0 
(i  = 1 , 2)  then  erfc  (h.)  is  required  where  I h£  I is  the  normal  distance  from  line(T)to  the  origin 
(See  (20)).  The  lines  (3)  and  J^4)  shown  in  the  diagrams  bound  the  angular  region  denoted  by  A.  In 
diagrams  [T),  [2],  (jQ,  A and  A coincide. 


If  I hi  ( I h j I or  I h2 1 ) is  sufficiently  small,  erfc  (h)  can  be  replaced  by  one  and  a call  to  the  erfc 
routine  avoided.  This  feature  appears  in  the  program  through  (1,12,13,16,20,34]  where  the  in- 
equality 

(44)  I hi  < 


is  tested.  We  have  if  (44)  holds 

(45)  j |erfe(h)  - l|  « ~ ihi  < otj  c «/2  , (e  defined  in  (13)) 
so  that  «3  is  taken  as 

(46)  s \/v  e/2. 


Box  (11,7 ) is  used  to  check  if  R is  sufficiently  large  for  the  computation  of  (19)  to  be  by- 
passed. The  choice  for  R A/3,  which  has  already  been  discussed  on  page  8,  is  made  so  that  with 
R > R, 

PIA(R,~~.~)1  < e. 

The  program  forP(A)  by  < 19)  is  displayed  »n{  11,12-23)  and  (II, 27],  with  |H,4)  showing  the 
computation  for  J0  which  denotes  t the  angle  of  A where  the  sign  agrees  with  the  sign  preceding 
P(A)  in  the  relations  given  for  P in  the  diagrams  of  Figure  5 (See  (18),  also). 

The  program  is  designed  to  recognize  and  avoid  a subtil  situation  that  can  occur  due  to  round* 
off  error  that  leads  to  a catastrophic  erroneous  result.  As  an  example  suppose  we  are  dealing  with  a 
polygon  where  one  of  the  exterior  angular  regions,  say  Ak . k =£  i.  N.  as  shown  by  the  solid  lines 
in  Figure  6,  subtends  an  angle  0 of  nearly  ? radians  with  sides  of  A at  large  perpendicular  distances 
from  the  origin,  so  tliat  P(A)  ~ i . 
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, * -• ' •*.}•  ■***.“  **•  -'.vu.  Sr  A 


Figure  6.  Shows  A Singular  Case  Situation 


Suppose,  however,  by  rounding  error  line  (7)  is  actually  given  by  the  computer  as  line  (5)  so 
that  the  angular  region  A subtends  an  angle  0.  near  (— it)  radians.  Thus  the  program  in  this  case 
would  yield  a value  ( ~P(X)) ; i.e.,  a small  value.  Moreover  it  would  be  negative  since  0 is  measured 
from  (3)  to  (I)  which  is  clockwise  rather  than  counterclockwise.  This  singular  case  situation 
and  others  that  can  occur  are  handled  in  the  program  through  boxes  (1,221,  (1,23),  (1,24),  and 
(1,291 . When  N > 3,  a singular  case  occurs  for  the  k,h  angular  region  of  H (A0  ( (0,  jr) ) when 


S s-^-sin<02 


0,)  ^ 0. 


If  this  is  the  case,  a second  inequality  is  tested,  namely. 


€ =-^-cos - 0 ,)  < 0, 


If  (47)  and  (48)  are  satisfied,  as  in  Figure  6,  we  set  P(Ak ) a Vi  erfc  (t),  where  t ® -hj  if  g,  < 0,  or 
t = h2  if  g,  > 0.  If  (47)  holds  and  (48)  does  not,  then  we  set  P(Ak)  = 0 since  I A0i  ^ I0"14.  A 
singular  situation  that  cannot  be  resolved  occurs  in  the  unlikely  case  that  (47)  holds,  (48)  docs  not, 
H > R,  ami  , g,  art-  negative.  When  ail  of  these  conditions  are  true,  Ak  may  contain  the  origin  so 
that  for  sufficiently  large  R (>  1 04 ),  P(Ak)  is  not  close  to  zero.  However  A0(*-  10"*4)  should 
always  be  in  (0,  rr)  for  a convex  polygon,  but  it  is  not  since  (47)  holds.  Hence  we  cannot  find,  with* 
in  the  single  precision  capabilities  of  the  CDC670O,  the  value  of  P(Ak),  because  the  value  of  A0 
cannot  be  resolved. 


In  the  next  section,  we  discuss  the  Gidcon-Gurland  method  (G  & G)  for  evaluating  P(A).  In 
their  report  and  published  paper  however  they  do  not  consider  the  programming  aspects  of  their 
method,  which  must  also  deal  with  the  singular  case  problem  just  mentioned. 

Extensive  checking  of  our  program  was  carried  out.  Comparisons  of  results  were  made  with  a 
program  of  the  G&G  method  that  we  developed.  Abo  comparisons  were  made  with  two  other 
independent  programs  for  computing  P(H)  for  the  special  case  of  triangles.  These  programs  also 
allowed  independent  checking  fer  convex  polygons  other  than  triangles,  since  a convex  polygon  can 
always  be  decomposed  into  a set  of  triangles. 
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Our  computing  program  is  designated  as  VALR16.  In  Appendix  C a Fortran  IV  listing  is  given 
of  a test  program  which  generates  coordinates  representing  the  vertices  of  a set  of  triangles  such  that 
all  phases  of  the  VALR16  program  are  tested  by  evaluating  the  probability  over  these  triangles. 

Some  numerical  results  are  also  giver,  there. 

5.  COMPARISON  WITH  GIDEON  AND  GURLAND  METHOD 

The  work  of  Gideon  and  Gurland  (G  & G)  (4),  [5]  gives  a set  of  relations  by  which  P(A)  can 
be  evaluated.  Their  unique  procedure  is  very  efficient  and  though  limited  to  5 decimal  digit  accuracy 
appears  to  be  one  of  the  best  of  the  methods  we  reviewed  in  the  literature  [ 1,  2,  3,  6,  7] . Essentially 
they  assume  the  angular  region  A has  been  rotated  as  shown  in  Figure  2 such  that  if  I0jl  < w/4, 

(i  = 1, 2)  then 

(49)  P[A(R,  0,  0j)l  =-^erfc(R/v/2)  [b0  + b,  R + b2  R2  )Qi  + (b3  R + b4  R2  )0 2 + (b5  R + b6  R2  )0f  )* 


The  coefficients  bj  were  determined  by  least  squares  for  each  of  IS  subintervals  in  R,  (0,  1/2), 

( 1/2,  .75 ) .... , (j/4,  (j  + 1 )/4) .....  13.75,  4) . In  order  to  evaluate  P(A),  they  need  to  use  (49) 
twice,  with  the  same  value  of  R,  once  fordj  and  once  for  02.  Because  the  use  of  (49)  is  constrained 
to  1 0 1 < ir/4,  G & G require  in  addition  to  (20)  the  relation 

(50)  P(A(R,  0,0)1  = ~erfc(^  sin  0)  crfc  (j~  cos  0)  - P[A(R,0,  y- 0)1  , ~ < 0 < */2. 

We  have  programmed  the  (G  & G)  method  and  found  the  average  computing  time  per  angular 
region  to  be  about  20%  longer  than  ours  at  the  6 decimal  digits  accuracy  level.  We  estimate  a 25% 
to  30%  difference  if  we  modified  our  method  for  5 instead  of  6 decimal  digit  accuracy. 

Although  it  takes  less  time  to  evaluate  the  righthand  side  of  (49)  twice,  without  erfc  (R/\/T), 
than  it  does  to  evaluate  the  recursive  procedure  given  by  ( 19),  our  method  lias  significantly  less 
calls  to  the  various  special  function  routines,  except  for  the  exponential.  In  particular  since  the 
minimax  approximation  for  erfc  (u)/z(u)  holds  for  u > 0.  ( 10  K ir/2),  we  do  not  need  (50). 
Moreover,  in  evaluating  the  number  of  erfc  functions  required  by  G & G it  is  recalled  that  we  need 
an  erfc  function  when  »/2  < 0 < 3t/4  or  when  3rr/4  <0  <rr.  fit  the  second  case  they  also  need 
one  erfc.  however  for  the  first  inequality  they  need  two.  Consequently,  for  each  A,  counting  the 
erfc  function  needed  in  (49)  once  and  using  (20)  and  (50)  it  is  easy  to  show  by  enumeration  of 
cases  (for  example,  in  (3j  of  Figure  5,  page  14,  they  could  neeu  5 while  we  would  need  none) 
that  their  method  takes  on  the  average  3Va  times  as  many  erfe  functions  as  ours.  In  addition,  they 
treat  0,  and  02  separately  while  we  treat  the  difference  02  - 0,  (except  for  the  functions  g, , 
h| , gj,  h3  which  are  expressed  as  algebraic  functions  of  the  coordinates  of  A).  Thus,  they  need 
two  separate  calls  to  tire  arctangent  routine  for  P(A)  whereas  we  require  one,  and  for  H a triangle 
they  need  5 arctangents  (taking  advantage  of  (37))  while  wc  need  only  2.  They  also  need  R which 
requires  a square  root  while  we  need  an  exponential.  The  average  number  of  calls  to  special  func- 
tions for  a convex  polygon  of  N sides  is  summarized  in  Table  1 . 

•We  note  a serious  omission  in  {$)  where  it  Is  not  explicitly  stated  that  (49)  only  holds  for  10 1 < */4. 
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Table  1.  Average  No.  of  Calls  to  Special  Function  Routines 
for  N-Angular  Regions 


Us 

G&G 

erfc 

N 

3.5N 

tan-1 

N-l 

2N-1 

square  root 

N+l 

2N+1 

exponential 

N 

0 

We  also  note  that  in  (4]  they  advocate  treating  N sided  polygons  by  decomposing  them  into 
sets  of  quadrilaterals  and  triangles.  In  the  case  of  N-convex  polygons,  this  would  mean  treating  2N-3 
angular  regions  for  N odd,  and  2N-4  for  N even,  whereas  we  would  only  require  N angular  regions  as 
explained  in  Section  3.  Also  in  the  case  of  an  arbitrary  polygon  it  will  be  more  efficient  in  general 
to  decompose  it  into  as  few  convex  polygons  as  possible  rather  than,  as  G & G propose,  into  tri- 
angles and  quadritalcrals. 


6.  COMMENTS  ON  DREZNER'S  METHOD 


In  a recent  paper  by  Drezner,  [3] , a method  is  given  for  computing  the  bivariate  normal  inte- 
gral. 4Km,  k,  p);  Le„ 

(51)  a (2rrv/T"'p-)“1  f f expl - \ ~ * ) dwdz. 

By  letting 

/ u “ (in  - w - p2)  , v » (k  - z)fy/2(\  -p3)’”. 

(52)  { _ __ 

(M  « mfs/Xr-  , R * kV2(l  - ph  . 


he  obtains 

(53)  <lHm.kj>)  = .'sffi-.r- j J e‘uJe'yJ  f(u,v)duuv, 

where 

(54)  Ru,v)  = cxp|M(2u  - M)  + R(2v  - K)  ♦ 2p(u  - MHv  - K)l. 
Drezner  then  uses  Gaussian  integration,  when  nt  < 0,  k < 0,  p < 0,  so  tltat 

JX oT  i i 

(55)  4Km.k,p)  a X— f-  v £ AAfto.v), 

» p*l  i»|  • > 11 


18 


where  the  weights  A(  and  abscissae  u(  (or  v.)  are  given  in  [8] . He  makes  the  significant  observation 
that  if 

(56)  m < 0,  k < 0,  p < 0, 

then  f(u,v)  < 1 and  he  can  use  (55)  directly  to  evaluate  <l>  within  a given  error  for  relatively  sma!' 
values  of  J.  For  example,  the  maximum  observed  error  for  J = 5,  is  reported  to  be  5,5(-7).  He  also 
takes  advantage  of  the  fact  that  if  the  argument  of  the  exponent  in  (54)  is  sufficiently  less  than 
zero,  f can  be  replaced  by  zero.  For  J = 5,  his  cutoff  value  is  stated  as  - 1 2. 

In  cases  where  one  or  more  of  the  three  inequalities  in  (56)  does  not  hold,  one  or  two  erfc 
functions  must  also  be  computed.  In  case  mkp>0,  then  two  sums  such  as  appear  in  (55)  are  needed 
in  addition  to  possibly  one  or  two  erfc  functions.  The  necessary  relations  are  all  given  in  (3) . Two 
typographical  errors  are  noted  there.  In  (10)  l/>/2ff  should  replace  \l2v  and  in  (12) 

„ 1 - Sgn(m)  Sgn(k) 

(57)  8mk  = 4 , 

where  the  minus  sign  replaces  an  incorrect  plus  sign. 

nearly  ( 1 ) with  p,  * 0,  o,  ® 1,  i =>  l,  2,  reduces  to  (51)  where  the  angular  region  A has  a 
right  angle  at  (m,k).  Applying  the  transformation  in  (2)  which  reduces  to 

(58)  x a (w  — pzj/v/l  - pi  . y a z, 

the  90°  angular  region  A in  tire  w - z plane  is  transformed  into  an  angular  region  A m the  x - y 
plane  with  vertex  at  (x0,  y0),  where 

x0  » (in  - pk)A/ 1 “pa  , y0  a k . 

with  a subtended  angle  0O . 

0O  » ta«“Vl  ~3"/(“P)l  . 

where  0p  is  measured  counterclockwise  from  tire  negative  side  of  the  line  y * k.  lire  angular  region 
A therefore  is  always  below  tire  line  y s k. 

In  particular,  given  a set  of  values  (rn.k.p)  there  exists  a corresponding  angular  region  in  the 
x - y plane  specified  by  R.  0, , 0j(  Sec  Figures  1 and  2).  The  connection  between  these  sets  of 
variables  can  be  shown  to  be 
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(59) 


R = l(m2  - 2pnik  ■>  k2 )/( 1 - p2)]54 

01  = tan-1  (k-v/l  - p2/(pk  - m)],  B~  ~ ian‘1[-r.'V»  ~ P2/(pm  - k)]  , 

gt  = (pk  - m)A/2(l  - pi-  g2  = (pm  - k)/^/ 2(1  - p2) 

hj  = kA/2,  h2  = -m/y/T , 


or 


m 


= -v/T“  sine2  = -v^h2,  k = v/I—rsingj  = s/lhx  (see  (40)), 


(60) 


P = -p(SiS2  + h^j)  = -cos (02  - 0j), 
0 - pT  - siu(02  - 6 x)  > 0. 


We  have  programmed  the  Dresner  procedure  and  compared  it  to  our  method.  A Fortran  IV 
listing  is  given  in  Appendix  B.  We  did  not  expect  it  to  be  as  efficient  because  of  the  large  number  of 
exponentials  required.  For  J - 5 25  exponentials  are  required  when  mkp  < 0,  and  50  are  needed 
when  mkp  > 0.  However  neither  method  suffers  in  comparison  to  the  other  in  computing  addi- 
tional erfc  functions  (or  livalently  normal  probability  integrals  in  one  dimension)  since  it  can  be 
shown  both  require  the  same  number  (none,  one,  or  two)  in  any  particular  case. 


Timingruns  -..r  the  two  programs  showed  that  the  Drezner  method  is  4 times  slower  on  the 
average  than  ours  for  6 decimal  digit  accuracy  and  8 times  as  slow  for  9 decimal  digits  of  accuracy. 

7/e  also  note  that  Drezner’s  procedure  is  incomplete  for  programming  because  he  does  not 
state  how  to  treat  tbc  cases  p = 1 , p = - 1 . These  values  can  occur  through  numerical  rounding  and 
must  be  dealt  with  before  a working  program  can  h ; obtained.  This  problem  is  resolved  by  noting 
that  if  p - 1 - e,  e > 0,  then 


(61)  lim  4>(m,k,  1 - e)  = ’Aerfcf-T l\/7), 
e -*■  0 


where  T = minimum  of  m and  k; 
ifp  =-l  +e,  e > 0,  then 


(62)  lim  4>(m,k,-l  + e)= 

e -*■  0 


'/2[erfc(-k Is/T)  - erfc(m/v'2)l , ifk>-m 


0 otherwise. 


These  formulas  are  easily  seen  to  be  true  by  noting  that  the  line  w = m transforms  by  (58)  to 
the  line,  call  it  L, 
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7 *•=! 


(63) 


x-y/l  - p2  + yp  = m. 


The  line  L dearly  has  the  property  that,  whether  m is  positive  or  negative,  it  is  tangent  to  the  drde 
(D) 

(64)  x2  + y2  - m2. 

The  following  additional  facts  an  easily  proved  from  (63)  and  will  help  the  reader  in  following 
Figures  7 and  8 and  similar  situations:  ( 1)  The  slope  of  line  L,  or  dy/dx,  is  + y/1  - p2/(-p),  and 
hence  the  slope  has  the  opposite  sign  to  that  of  p,  so  that  p must  be  negative  in  Figure  7 and  positive 
in  Figure  8;  (2)  The  x-intercept  of  L is  m//l ! - p2,  so  that  m is  positive  in  Figure  7 and  negative 
in  Figure  8,  m having  the  same  sign  as  the  x-intercept  of  L. 

Hence  if  p = -1  + e,  e > 0 and  small,  k > -m,  we  have  the  situation  shown  in  Figure  7,  Now 
as  e -*>  0,  the  point  (C)  approaches  + <»  along  y = k,  and  L approaches  tangency 


to  the  circle  (D)  at  (O.-  m)  (but  note  that  0 < k < m in  this  case  since  the  x -intercept  of  line  L is 
positive).  Consequently,  in  the  limit  as  e 0,  <Km,k,-  l)  is  given  by  (62).  Similar  heuristic  argu- 
ments. which  can  be  made  rigorous,  can  be  given  for  any  other  situation  with  p -»  1 as  well  as 
p -*  “1.  For  example,  with  p 3 1 - e,  e > 0,  m < k as  in  Figure  8. 


In  this  figure,  p = 1 - e,  m < k < 0 (x-intercept  negative),  and  we  observe  that  as  e -*•  0,  the 
shaded  area  in  Figure  8 approaches  the  region  below  and  including  the  line  y = m(<  k < 0)  as 
required  by  (61),  or  the  limit  is  ( 1/2)  erfc(-m 


7.  SOME  NUMERICAL  RESULTS 


in  this  section  we  give  the  numerical  results  for  P(H, ) and  P(H2),  using  our  program  VALR-16, 
where  H(  is  a 6 sided  convex  polygon  containing  the  origin  and  H2  is  an  8 sided  convex  polygon 
not  containing  the  origin.  The  x,  y columns  of  Table  2 below  contain  the  x,  y coordinates  of  the 
vertices;  the  three  columns  that  follow  list  the  values  of  P(Ak ) for  each  angular  region  (See 
Figure  3)  for  6j  a 2.5(-4),  e2  as  2.6(-7),  e3  a 2.9(-10).  The  last  row  headed  P(H)  contains  the 
value  of  P(H)  for  t't , e2 , e3 . Al!  the  P values  have  been  truncated  from  14  digit  CDC  6700  output. 


Table  2 


k 

— 

X 

y 

P(Ak),e, 

P(Ak),  e2 

P(Ak),  e3 

1 

o.s 

-2.0 

.911227 

.911064879 

.91106477067 

2 

2.0 

0.0 

.046858 

,046998988 

.04699911886 

3 

0.5 

2.0 

.052500 

.052666886 

.05266699792 

4 

-0.5 

1.5 

.059487 

.059482771 

.05948276788 

5 

-1.5 

0.0 

.042640 

.042515227 

.04251511748 

6 

-1.0 

-1.5 

.017780 

.017747368 

.01774728061 

PUD  = 

P(H,)- 

.727S21 

.727148375 

.72714804914 

1 

1.5 

-1.5 

.851151 

.850975856 

.85097578896 

n 

2.0 

-0.75 

.018552 

.018726585 

.01872664573 

3 

1.75 

1.75 

.038192 

.038305815 

.03830565474 

4 

1.25 

1.^5 

.064039 

.064042498 

.06404250011 

5 

0.50 

1.50 

.486789 

.486841686 

.48684172637 

6 

0.25 

0.25 

.042253 

.042256194 

.04225618944 

7 

0.25 

-0.25 

.026270 

.026273494 

.02627348809 

8 

0.50 

-1,25 

.116865 

.116775266 

.11677520430 

PUD  « 

P(Hj)- 

.291919 

.291304849 

.29130478878 

(See  (36)) 
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APPENDIX  A 

PROGRAM  PARAMETERS.  CHEBYSHEV  COEFFICIENTS  FOR 

erfc(x)/z(x),  x > 0 

« 

In  this  appendix  we  list  the  pertinent  constants  that  appear  in  the  program  for  three  levels  of 
accuracy  (3,6,9  decimal  digits),  and  an  additional  set  which  is  designed  to  yield  12  correct  decimal 
digits  for  the  probability  over  an  angular  region. 


6 

a«)  ^ 

€ 

al 

“2 

«3 

yE(1A/D 

4.50(-4) 

2.46 

2.54(-4) 

2.02(-7) 

1.22.(— 2) 

2.25(— 4) 

2.52(— 4) 

4.56(-7) 

3.5505 

2.57(— 7) 

2.080-13) 

1.23(-4) 

2.28( — 7) 

2.57(— 7) 

5.2K-10) 

4.382 

2.94(-10) 

2.72(— 1 9) 

1.35(— 6) 

2.6K-10) 

2.88(— 10) 

l.78(-13) 

5.1092 

l .00(— 1 3) 

3. 1 7(— 26) 

6.58(-9) 

8.90(— 14) 

2.50H3) 

e = S/y/n' 
C<6) 

KA/T 


a,  = ire* 


See  page  6. 

See  page  5. 

See  page  8. 

See  pages  7,  1 3. 


a2  = (9ire2),/3  See  pages  7, 13. 

ofj  = 6/2  See  pago  15. 

- E(R'A/2)  — ~erfc(R/N/2)  See  page  8. 

= 2.5c  (for  ® ) 


The  first  column  of  the  table  labeled  Acc.  (for  accuracy)  lists  0 .(§)>©*  (§)  rcferring.to 
3,  6. 9.  1 2 decimal  digits  of  accuracy,  respectively,  for  the  probability  over  an  angular  region.  Pages 
are  given  above  where  the  parameters  arc  defined  in  the  report. 


The  minimax  coefficients,  ak,  for  approximating  erfc(x)  on  C(5)(See  (12),  (15))  are  given  be- 
low for  four  accuracy  levels  as  indicated  in  the  tables  below  by  (X)  ,(§),©,  (tj)  . They  were 
computed  by  a double  precision  minimax  subroutine  utilizing  values  of  erfc(x)  accurate  to  18  sig- 
nificant digits  on  ('A,  C]  and  erf  (x)  accurate  to  25  digits  on  (0,V5) . 

For  0 (Average  time  per  angular  region  a 2.2  x 10‘4sec) 

^ = .8857775 18572895D  + 00  a,  « -.981 ! 5 1 9527780S0D  + 00 

aj  a ~ .35 3644980686977D  + 00 


a2  » .759305502082485D  + 00 
a4  - .695232092435 207D  - 01 


For  (b)  (Average  time  per  angular  region  = 4.6  x 1 O'4 sec) 


a0  = .8862264700 16632D  + 00 
a2  = .885348820003892D  + 00 
a4  = .42182 1 197 160099D  + 00 
a6  = .905057384 150449D  - 01 
ag  = .4308951 68984 138D  - 02 


For  (C)  (Average  time  per  angular  region  = 


a0  » .886226924931 465  D + 00 

a2  = .886223733 186722D  + 00 

a4  = .44285 1899328568D  + 00 

a6  = .14506004340301 2D  + 00 

a8  = .309199295521210D  - 01 

at0=  .324944543171 185D  -02 
a,  2=  . 105787574480633D  -03 
al4  = .4083355 17232165D  -06 

For  (§)  (Average  time  per  angular  region  = 


ao  • .886226925452593D  + 00 

a2  = .886226922786746D  + 00 

a4  = .443112868048919D  + 00 

h a .1476871363219380  + 00 

a8  » .3680328493508600  - 01 

aio=  .7 10292625734052D  - 02 
at2=  .98111 2629090333D  - 03# 
UI4“  .789960968802448D  - 04 
ai  6“  .2836466354093220  - 05 
al8=  .3176794970400060  - 07 
a2o B .452S34347337305D  - 10 


a,  = -.999950714561 036D  + 00 

a3  = -.66061 1239043357D  + 00 

a5  = -.222898055667208D  + 00 

a7  = -.2549061 1 1884287D  - 01 

a9  = -.323377239693247D  - 03 


.5  x 10"4  sec) 


a,  = -.999999899776252D  + 00 
a3  = -.6666266705 10907D  + 00 
as  = -.265638206366O25D  + 00 
a7  = -.7 14909837799 889D  - 01 

a,  = -.112323532148441D  - oi 

ai  i = -.704260243309096D  - 03 
al  3 - -.9718648641604610  - 05 


.1  x 10-4  sec) 


a,  « -.9999999999485970  + 00 
a3  » -.6666666 11 86666 ID  + 00 
as  n -.2666627290914110  + 00 
a7  = -.76136585S850292D  -01 
a9  « -.1671950968881830  - 01 
al  t°  -.278 1 7Q932906224D  -02 
an  = -.302S88640752108D  -03 
ais“  -.1686851817670460  - 04 
al7=  -.358314466908290D  -06 
ai 9“  -.I7S4406S1940430D  -08 


Average  time  per  angular  region  refers  to  the  average  computing  time  on  the  CDC-6700  to  obtain 
P(A). 
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APPENDIX  B. 

LISTING  OF  DREZNER  PROGRAM 

This  appendix  contains  a listing  of  the  program  for  computing  P(A)  or  ?(H)  by  Drezner’s  procedure. 
It  is  designed  to  use  J = 3,  5,  8 where  J is  defined  by  equation  (5)  in  [3] . Thus,  referring  to  Table 
1 in  [3] , P will  be  computed  correctly  to  at  least  3, 6.,  or  9 digits,  respectively  by  this  program. 


Call  line  to  Z.  Drezner  Subroutine 

CALL  DREZNR  (x,  y,  N,  Pk , IOP)  where 

x is  the  input  array  of  abscissas  of  the  vertices  of  the  polygon  , , , . 

I Vertices  must  be  listed  in 

y is  the  input  array  of  ordinates  of  the  vertices  of  the  polygon  ' 9ounterclock\ifise  order.  See  pp.  9, 10. 

N is  the  number  of  sides  of  the  polygon.* 

Pk  is  the  location  of  the  answer  as  computed  by  the  Drezntr  method 

IOP  « 1 specifies  the  Drezner  subroutine  to  use  a 

table  of  J = 3 weights  in  computing  Pk  (Sec  (55)). 

IOP  2 specifies  the  Drezner  routine  to  use  a table 
of  J = 5 weights  in  computing  Pk. 

IOP  = 3 specifies  a table  of  J = 8 weights  in 
computing  Pk . 


*N  = 1 for  an  angular  region  A with  3 points  given  in  counterclockwise  order  with  first  point  at 
vertex  of  A,  (See  page  12).  Note  0 < AO  < 2v  for  N»  l.butO  < AO  < srforN  > 3. 
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. V vrjiB:  - 


SUBROUTINE  OR£ZNR(  X, Y,N,ANS,  IOP  ) 

OIHCNSXON  X(l)»Y(l)»U(2)fV(2)»G(2)»H(2) 

OIMEN SION  AM(51),AK(5i),RH0(51) 

REAL  L 

OATA  RT2  / 1.4142  13562  373  / 

NMl=N-i 

K=1 

ANS=0, 

NBAR=N 

IF  ( N.EQ.i  ) NBAR=3 
U(2)=X(NBAR)-X(1! 

V(2)=Y(NBAR)-Y(1) 

KPi=K+l 

U(l)=X(KPll-X(K) 

V f 1 Is YCKP1I-Y (K1 

IF  ( N.GT.l  I GO  TO  3141 

SGN=1. 

SNsV(2l*U(l)-U(2)*V(l) 

XF  ( SN.GE.O.  ) GO  TO  3141 
SGN=-  1. 

T 1*  U ( 11 
U ( 1 )“U ( 2) 

U(2  ) = T1 
Tl= V ( ?> 

V(2)=  V( 1) 

V ( 1 I=T1 
3141  CONTINUE 

BGQ1*SQRT<  2. ♦<U< 1>*UC1» ♦Vll ) *V Sill! 

8G02*SQRT ( 2.*  (U(2>*U  (21  ♦V(2)  *V  (211) 

3151  CONTINUE 
1 = 0. 

B = .5* (X (Kl *X ( K»  >Y  <K»  *Y (K I > 

G(t>*Ofi)*X(K>4V(l>*Y(K) 

G(2)=U(2l*X(KltY(2)*Y(K) 

M(l)«-V(K*»UI1I#X(K»*VC1I 

H(2I=-Y(K>*U(2)4X(KI*V(2) 

6(1  >=€  (D/BGD  1 
G(2 )=C(2)/BG02 
M(l)=Mi)/BG01 
H(2I  = M 2I/BG02 
AH(KJ=-RT2*H(2i 
AK(K) *RT2*H(1 ) 

IF  ( F.NE.0.  I GO  TO  3131 

RH0(KI*»(2. MU(2I*U(1I+V  (2>*V  (1 )) )/ (BG01*6G02) 
GO  TO  3191 
3101  CONTINUE 
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RH0<KI=-<Gfi**6(2H-HUI*M2H/B 


S ' 

ft 


P 


3191  CONTINUE 

IF  ( K.LT.NMi  ) GO  TO  3631 
IF  ( K.EQ.NMi  I GO  TO  3661 
CALL  FLAN  ( AM ( K) , AK ( K) , RHO (K ) ♦ ANSI » IOP  ) 

ANS=A  NS+ANSl 

IF  < N*NEo 1 > RETURN 

IF  ( SGN.EQ.l.  ) RETURN 

ANS=1.-ANS 

RETURN 

3631  CONTINUE 
K=K+1 
KPi-K+i 

IF  ( K.NE.2  ) GO  TO  3651 
KHlaK-1 

CALL  FLAN  ( AM ( KM  1) , AK (KH 1)  , R FO  !KH1 ) , ANS  1,  IOP  ) 
AN S= ANSI 
U(2»=X<KP1I-X  <K> 

Vt2l=Y«KPl)-Y(K> 

BG02=SQRT(  2.  ♦ <U  < 2I*U(2J  +\li  2)  *V  C2M> 

GO  TO  3151 
3651  CONTINUE 
Uti )=U<2) 

V < 1 )=  V<  2) 

U«2l=VtKPll»X(K» 

V«2»=Y(KP1)-Y (K» 

9601=  EGO? 

BGO?=  SQRT < 2«*<U(Z>*U(2l*\M2l#\M2IH 
GO  TO  3671 
3661  CONTINUE 
K=N 

U*II»X<N»-X<1S 

V<n«YtN)-Y(lt 

9G0i=SQRT(  2*  *CU<1I*U(1MV<1>  *V  <ll)  1 
3671  CONTINUE 
KH1*K *1 

CALL  FLAN  t AM(KM1>  » AK  (KN1 ) «RN0CKH1 ) t ANS  it  I OP  ) 

ANS*ANS-ANS1 

GO  TO  3151 

ENO 
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SUBROUTINE  PLAN  ( H,AK,R,ANS,  IOP  ) 

DIMENSION  EPS 3 (11) 

OATA  C EPS3<I ) f I=l»3  ) / 2.E- 5»  2.E-7*  2.E-1Q  / 
OATA  RT 2/1*4142  13562  373  t 
0M=l*-EFS3tI0P> 

ANS=0. 

IF  ( R.LE.-OM  ) GO  TO  3171 
IF  ( (H*AK*R) .GT.O.  ) GO  TO  3155 
IF  ( H.GT.Q,  ) GO  TO  2031 
IF  ( AK.GT.Q*  ) GO  TO  2021 
IF  ( F*GT*0»  ) GO  TO  2011 
ANS-8FHI (Ht AK*R* IOP  ) 

GO  TO  3161 
2011  CONTINUE 

IF  ( AK.NE.O.  ) GO  TO  2061 
GO  TO  2023 
20  21  CONTINUE 

IF  ( F.LT.O*  ) GO  TO  2041 
20  23  CONTINUE 

ANS=EC9(H,AK,R,I0P  ) 

GO  TO  3161 
20  31  CONTINUE 

IF  ( AK'EQ'O.  ) GO  TO  2051 
20  35  CONTINUE 

IF  ( AK.LT.O.  > GO  TO  2061 
2041  CONTINUE 

ANS*EC7 (HiAKf R»IOP  ) 

GO  TO  3161 
2051  CONTINUE 

IF  ( F.CT.Q,  ) GO  TO  2061 
GO  TO  2041 
2061  CONTINUE 

ANS*EC8(H,AK*R,I0P  ) 

GO  TO  3161 
3155  CONTINUE 

ANS*E  Cl 1(H*AK#R« IOP  ) 

3161  CONTINUE 
RETURN 

3171  CONTINUE 

IF  ( AK.LE. (*NN£PS3 (IOP) 1 1 GO  TO  3161 

T1*»A  N/RT2 

T2*H/RT2 

ANS». 5*<ERFC(O*Tl)-ERFG<0tT2l ) 

GO  TO  3161 
ENO 
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FUNCTION  EOT  (H*  AKtR* IOP  ) 

OAT  A RT2/1*4142  13562  373  / 

T*-H/PT2 
Tl*-A  K/RT2 

EQ7=BPNIC-H»-AK*R*I0P  »4*.5*<ERFC«0»TS+ERFC<0*Tin-l. 

RETURN 

ENO 


FUNCTION  EQ6  <H,AK»R,IOP  ) 

DATA  RT2/1.4142  13562  373  / 

T«-AK/RT2 

EQ«*-EPHI<-H*AK*-Rf  IOP  > ♦•5*ERFC«0,T) 

RETURN 

ENO 


FUNCTION  EQ9  <H,AK.R,IOP  ) 

DATA  RT2/1*  4142  13562  373  / 

T*-H/RT2 

EQ9«-BPH<M,-AK»-R,I0P  ) ♦•5*ERFC«G.T> 

RETURN 

ENO 


31 


FUNCTION  EQ11(H, AK,R»lOP  > 

OIMENSION  EPS3(11) 

OATA  C EPS3(I  i ,1*1,  3 > / 2.E-5, 2.E-7, 2.E-10  ✓ 
DATA  FT 2/1* 4142  13562  373  / 

91  FORMAT  ( 1H0,  3E22»15  ) 

ON*1,-EPS3(IOP> 

IF  ( fi.LT.OM  ) GO  TO  2001 
T = H 

IF  C AK.LE.H  ) T*AK 

Ti=-T/RT2 

E011*,5*ERFC(0,T1) 

GO  TO  1991 
1991  CONTINUE 
RETURN 

2001  CONTINUE 

CST  *SQRT (H*M»2, ♦R*H*AK^AK*AK  ) 

T1=R*H-AK 

Cl=l, 

T2*SI GN  (Cl, Hi 
Tl*(Ti«T2)  /CST 
T4=l. 

T3=H*AK 

T5*SI CN  (T4»T3I 

TOELXlf  T5)  •* 25 

T3=P* AK-H 

Cl»l. 

T2=SIGN (Cl, AK) 

T3=(T3*T2>  /CST 
IF  ( P.GT.O.  » GO  TO  2031 
IF  < T1.GT.0.  ) GO  TO  2023 
T4*0PHI(M,O.,T1,IOP  ) 

GO  TO  2051 
20  23  CONTINUE 

T4*EQ«CM,Q.,Ti*IOP  » 

GO  TO  2051 
2031  CONTINUE 

IF  t T1.LT.3.  > GO  TO  2041 
T4*EQe<MtO**Tl,IOP  J 
GO  TO  2051 
2041  CONTINUE 

T4#E07IH,0,,T1»IOP  I 
2051  CONTINUE 

IF  I AK.GT.O.  I GO  TO  3031 
IF  ( T3.GT.0.  I GO  TO  3023 
f 6s0PNI  (AK,G« ,T3, I OP  > 

GO  TO  3351 
30  23  CONTINUE 

T6=EO«(AK,0. »T3,I0P  } 
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3031 


3041 

3051 


1 

2 

3 

4 

1 

1 

3 

4 

1 

2 

3 

4 

1 

2 

3 

4 


GO  TO  3051 
CONTI NUE 

IF  { T3.LT.0.  > GO  TO  3041 
T6»EOe(AK,0.»T3,IOP  ) 

GO  TO  3051 
CONTINUE 

T6*EQ7(AK,0.,T3,I0P  ) 

CONTI NUE 

EQll*T4*T6-T0Et 

RETURN 

ENO 


FUNCTION  dPHl  ( H,AK,R  , IOP  ) 

DIMENSION  A<21),X121),U0(6)  ,IHI (6) 

DIMENSION  EPS1 (11  ) 

DIMENSION  EPS3 < 11 ) 

DATA  ( Ad)  .1=1,8  ) i 
4.4602  9770  4 66658E-1, 

4.3728  88798  77644E-2, 

3.9233  lc 666  523996-1, 

3.3246  66035  134396-2, 

OATA  < X (II,  1=1,6  ) / 

1.9055  41497  96192E-1, 

1.7997  76578  415  73E*0, 

4.8261  39660  46201E-1, 

1.7797  29418  52026E*Q, 

OATA  < A(I), 1=9,16  ) / 

1.3410  91884  533606-1,  2.6833  07544  7264QE-1, 
2.7595  33979  864226-1,  1.5  744  82826  1879QE-1, 
4.4814  U991  746256-2,  5.3679  35756  ( 2526E-3, 
2.0206  36491  324376-4,  1.1925  96926  595326-6  / 
OATA  ( X (. ) *1=9, 16  ) / 

5.2976  64393  185146-2*  2.6739  83721  677676-1, 
6.1630  28841  62402  E-l,  1.0642  46312  116236*0, 
1,5608  55862  270066*0,  2.1639  21153  395866*0, 
2.8631  33883  708086*9,  3.6660  07162  724406*0  / 
DATA  ( EPSKI.'  ,I»1,3  ) / -8  .,-12., -20.  / 

OATA  PI  / 3.1415  92C53  58979  / 


3.9646  62669  98335E-1, 
2*4840  61520  28443E-1, 
2.1141  81930  760576-1, 
a. 2485  33445  15628E-4  / 

8.  4825  lo675  445772-1, 
1.0024  21519  682166-1, 
1.0609  49621  525726*9, 
2.6697  60356  36766E*0  / 


OATA  t LLUil ) ,1  = 1 ,3  ) / 1,4,5  / 

OATA  < LHI  ( 1)  , 1=  1, 3 ) / 3,6,16  t 
OATA  RT2  / 1.4142  13562  373  / 

OATA  t EP33U),  1=1,3  ) / 2.6-5,2.6-7,2.6-10  / 
0M*1.-EPS3(I0P) 

II 0*110 1 IOP) 

IHI=LHICIUP) 

EPS=EPSlUOP) 


RSQ*R*R 

IF  I RSO.Lt . 1.  ) GC  TO  2991 
T3*  l, 

CST=0 . 

GO  TO  3001 
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•v  ' ;v-V* \ 


2991  CONTINUE 

T3*SQRT ( l.-RSQ) 

CST=RT2*T3 
3001  CONTINUE 
8PHl=0. 

IF  < R.LE.-Ort  ) 60  TO  3011 

IF  ( ft.LT.01  I GO  TO  3331 
T=H 

IF  < AK.LE.H  J T-AK 

T1=-T/RT2 
8PHI=»5*ERF3 (0»T1) 

GO  TO  3371 
3011  CONTINUE 

IF  ( AK «LE. C~H*EPS3 ( IOPI  » ) GO  TO  3371 

T1=-AK#RT2 

T2=H/RT2 

ANS=.5*{ERF3<0,Tll-ERFC<0,T2n 
GO  TO  3371 
3331  CONTINUE 
H1*H/CST 
AK1*  AK/CST 
SUN=0. 

00  3361  1=1 LO# INI 

SUN 1*0. 

00  3351  4*11.0# INI 

T 1«H1*  12  .*XC  XI  *»H1  J*AK1*<2.*X  t J}  - AK1) 

1 ♦2«*R*{XUI-Hl»MXUJ-AKll 
IF  < T1.LT.EPS  I GO  TO  3351 
5UHl*SUMl*£XP<Til*A<0) 

3351  CONTINUE 

SUH=SUH*A<I»*SUH1 
3361  CONTINUE 

aP«I*ISUH*T3l/PI 
3371  CONTINUE 
RETURN 
ENO 
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APPENDIX  C. 

LISTING  OF  TEST  PROGRAM  WITH  SOME 
NUMERICAL  RESULTS 

Tile  test  program  listed  in  this  appendix  is  designed  to  “see”  all  paths  of  the  VALR1 6 sub- 
routine, the  basic  routine  of  this  report.  The  test  program  treats  three  different  sets  of  triangles  for 
each  e.  i.t\,  e, , e2 , e3.  A total  of  35 1 triangles  are  treated.  Our  subroutine  VALR16  is  used 
to  obtain  the  probability  P(H)  over  each  triangle  and  the  result  is  compared  with  the  result  obtained 
by  the  routine  based  on  Dreznefs  method.  The  numerical  results  below  state  the  case  number, 

(x,  y)  vertices,  VALRio  result,  Drezner  result,  and  absolute  value  of  the  difference,  corresponding 
to  that  case  for  which  the  absolute  value  of  the  difference  in  P(H)  for  the  two  methods  was  a maxi- 
mum for  each  set  and  for  each  e.  Thus  there  are  nine  cases  given  below. 


P(H)  and  I A PI 


e,  = 2.54(-4)* 

3 

2.0000  00000  0000 
1.0000  00000  0000 
3.0000  00000  0000 

1.0000  00000  0000 

0.0000  00000  0000 

1 .0000  00000  0000 

.01116  23895  4828 
.01144  55124  4546 
2.83(-4) 

116 

3.0000  00000  0000 
0.0000  00000  0000 
3.000.:  00000  0000 

1.5000  00000  0000 
-0.0006  06881  7000 
0,0000  00000  0000 

.07276  76379  5214 
.07312  88147  1695 
3.61(-4) 

76 

.11048  34376  7180 
-1.8895  16562  3282 
-1.8395  16562  3282 

.11048  34376  7180 
.1  1048  34376  7180 
-.88951  65623  2820 

.07464  96837  3470 
.07443  77215  8773 
2.12(~4) 

e2  = 2,57(~7)* 

15 

0.0000  00000  0000 

1 .0000  00000  0000 

0.0000  00000  0000 

-2.0000  00000  0000 

0.0000  00000  0000 

1 .0000  00000  0000 

.17865  07387  5631 
.17864  99501  5890 
7.89(-7) 

90 

3.0000  00000  0000 

3.0000  00000  0000 

6.0000  00000  0000 

3.0000  00000  0000 
0.0000  00000  0000 
1.5000  00000  0000 

.00059  65636  6379 
.00059  75925  2985 
9.29(~7) 

79 

.01109  91882  5761 
-1.9889  00811  7424 
.01109  91882  5761 

.01109  91882  5761 
.01  109  91382  576! 
-.98890  081  17  4239 

.11700  59368  5515 
.11200  54478  9902 
4.89t-7) 

e,  = 2.94(-10)* 

96 

4.0000  00000  0000 

2.0000  00000  0000 
0.0000  00000  0000 

0.0000  00000  0000 
2.0000  00000  0000 
7.8256  90500  (-10) 

.12383  33015  1674 
.12383  33012  5254 
2.64(-10) 

65 

1.9999  80000  0000 
-3.0000  00000  0000 
3.0000  00000  0000 

4.5000  00000  0000 
0.0000  00000  0000 
0.0000  00000  0000 

.47645  25380  3718 
.47645  25375  5211 
4.85(-10) 

: 94 

.00116  10499  1180 
-1.9988  38950  0882 
2.0011  61049  9118 

.00116  10499  1180 
-1.9988  38949  3056 
-1.9988  38950  0882 

.22803  35273  2106 
.22803  35268  0156 
5.19(-10) 

PROGRAM  DRET  ( OUTPUT  ) 

COMMON  IOP 

DIMENSION  X ( 201)  ,Y ( 201) ,X1(2Q1> »Y1 (201) 

OINENSICN  EPS1(4)»X3(3),Y3(3) 

OIMENSI  CN  APH2K3)  ,APH31<  3) 

DIMENSION  IRAY (21) 

OIMENSI  CN  OEt  1(  o ) » 3 1(3) , ALPHA  ( 3) 

C SET  A FOR  PJ  )/  HULL  OEGK 

DATA  < Xd),I=t,48  ) f 

1 l.,3.,2«,  1.,  2.  , “l*  , 1«,2.,**1«,  1.  ,0.  » 0*  t 

2 1«,0.,0«,  1.  , “1. ,2. , l.,0.,P.,  1 • , 0.  , 2 «, 

3 1*,Q.,3«,  l«)2t|  2«|  i.,2.,2.,  l.,3«,Q«, 

% Q«,2*,l.,  0»,i»,*i*»  0«f  “1  . ,1  * , 0., 1* , 2.  / 

DATA  ( Y(I),  1=1,  48  ) ( 

1 0«,1«,1«,  0.,1«,1«,  Q*,1*9~1«,  0*,2.,1«, 

2 0*, l.,-2.,  0*t 1«, -2.,  0»»  •i.,»2.»  0.»-l*» *1.  » 

3 0 » ,-2* , I*;  0 • , »2.»  -i.  , 0.  , -1., 2.,  0 .,  *1.  , 2.  , 

4 0«»i*»l.»  0»,1.,1«,  D.,-2.,1.,  0.,-l.»*l»  f 

C SET  B FOR  PJ  SV  HULL  3EGK 

DATA  ( X(  I)  ,1*49,90  ) / 

1 1«,«5,2«,  1. , 1. ,2.,  1 ., 1.5,2.,  l.,.5,'l., 

2 li,li,*l«,  1. , .66 666 , ~1* , 1.  , 2.  , 1.  33333 , 1*,2.,1«, 

3 1.  ,2. ,.5,  l«,”l«i,“l»,  1* , *1*  ,1*  , l«,lt5,*  li  , 

4 l.,l.,2«,  i.,2.,1.  / 

DATA  ( Yd), 1=49, 90  ) / 

1 0., *1.5,0.,  0 « ,* 1. , 0.  , Q • ,-• 5, 0 • , 0.,«5,0«, 

2 0 • , 1«,  0 • , 0.,  1 .5,  0 •,  tJ.,0.,.5,  Q.,0.,1«, 

3 0 . ? 0 . , 1.  5,  0 • , 0 • , **  1* , C. , 0 •,  -1.  , 0«,1«,Q.  , 

4 0*  ,*  2*  ,1  • , 0 • , .5 , 1.  t 

OATA  ( X(I),  1=91,102  ) / 

i .5, 1.5, 0«,  2. , 0.  , 4«,  2.,Q.,4«,  .5, 1.5,0.  / 

OATA  ( Y ( I)  , 1=91,102  ) / 

1 .5,1. ,2.,  2* , 1 * , 0 • , 2«,i.,0«,  *5,1*,  2*  t 

OATA  ( XC I),  1=103,  117  ) / 

1 1 • , 0 » ,0*  , 1. ,0  « ,Q . , 1 » ,0  * , 0. , 1*,Q.,1., 

2 l.,l.,G.  / 

OATA  ( Yd), 1=103,117  ) / 

1 0.,i«,l.,  Q.»l.»*.5,  D.,«5,l., 

2 0.,1.,*.5,  0*i*9, 1*  / 

OATA  i EPSKD  , 1=1,  3 ) t 

1 2.5362  66450E-4,  2.5714  45247; -7,  2.  9434  48712E-1Q/ 

OATA  ( DFLld),  1 = 1,3  ) t 
1 4.49542E-4,  4.55  777E-7,  5.217127E-10  / 

OATA  < Bl(I), 1=1, 3 ) / 2.46,3.5505,4.332  / 

OATA  ( APH3KI)  ,1=1,3  ) / 

1 .22477l£*3,  » 2278  8 85E-6,  • 26Q85635E-9  / 

OATA  ( APH21d)  ,1=1, 3 ) / 

1 .12206  59E-1,. 12319  198F-3,  .13480  369E-5  / 
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90  FORMAT  ( 1H0,60K,I10  ) 

91  FORMAT  < 1H-,8£16.8  > 

92  FORMAT  < 1H-  ) 

93  FORMAT  ( 1H0,I10  > 

97  FORMAT  ( 1H1  ) 

95  FORMAT  ( 1H  ,44X,2F12.4  ) 

N*3 

PRINT  97 

DO  3021  1*1,151 

xm>*x{i) 

Yl(I)sY<I) 

3021  CONTINUF 

00  3071  13*1,3 

PRINT  97 
PRINT  90,13 
IFNT*  0 

OO  3061  15*1,3 

PRINT  92 
TMAX*  0* 

M*0 

OO  3025  1*1,151 

X(I)aXlCI) 

Y(X)*Yi(I) 

3025  CONTINUE 

APM3*  APH31( 13) 

Y(92)*1.5-3,*A»H3 

Y(95)=3.*AFM3 

Y (98)  * -3»* APH3 
Y(101)*l • 5*3*  *49H5 
Y(lQ4l*.9*APM3 

Y (1  05 ) *-  • 9*APH3 
Y ( 107) sY (104) 

Y (111) =Y (105) 

Y (113) *Y (104 ) 
Y(ll7)*Ya05) 

OO  3027  t*t,l>l 
Y1 (I) ®Y( I) 
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3027  CONTINUE 

IF  { 1 5*EGU1  ) GO 
IF  ( I5.EQ.3  ) GO 
DO  3031  1*1 ,151 

X(I  ) = 3**X1(I ) 
Y(I)=3.*Y1<I> 

3031  CONTINUF 
3035  CONTINUE 

CO  3051  1=1,115,3 

3037  CONTINUE 
IP2=I*2 

DO  3049  11*1,3 

M=M*i 

IF  ( IPNT.NE.O  » I 
3041  CONTINUE 

IF  ( Ii.EG.l  ) GO 
T 1*X ( I) 

T2*Y(I> 

X(  I)  *X  (I  Fl) 

Y(I)  = Y(I*1) 
X<m>*XIF2> 

Y < I ♦ 1 ) *Y (IP2) 
X(IP2)*T1 
Y(IP2)*T2 
3047  CONTINUE 

IF  ( I5.NE.3  ) GO 
IF  ( Il*r,T,i  ) GO 
T1*APH21<I3> 

T 2*SQPT(T1|  -1*  £•  12 

$X*X(I)-T2 

SY*Y(I)-T? 

00  3043  17=1,1*2 

X(I7)*X(X7)-SX 
Y(I7)  = Y(I7)-SY 


GO  TO  3035 
GO  TO  3 035 


PRINT  93,  M 
TO  3047 


GO  TO  3045 
GO  TO  3045 


3043  CONTINUE 
3045  CONTINUE 

3048  CONTINUE 

T1»(X(I)-X(I*1)  ) *<V(I)-Y(I+2)) 

T2=(  X ( I)  -X{I+2)I  *(Y  (I)  *Y  (I+i) ) 

IF  < T1.GE.T2  ) GO  TO  304 
T i=X(  1+  2) 

X (I  +2)  s > (I  +1 ) 

X < I+i ) *T1 
T2=Y  ( I+?) 

Y(I+2)sY(I+i> 

Y(  1 + 1)  *T2 
304  CONTINUE 

IF  < IFNT.NE.O  ) 

t PRINT  95,(X(J),f  <J>  ,JsI,IP2  ) 

I0P1=I3 

CALL  ALR16(  X (I ) , Y (I) ,N, ANSI »I0P1  ) 

ANSsANSl 

IOP*I3 

CALL  DRE2NR  ( X{  I)  , Y<I>  ,3  ,ANS3,  13  ) 
OEL® A8S ( A NS- AN33) 

IF  ( OEl.LE.TM4X  ) GO  TO  3049 

MSAtf=  N 

SAVDEL=GEL 

TNAXs  DEL 

X3( 1) «X( I) 

X3(2)sX<I+i> 

X3<3)  = X(I+2> 

Y3(  1>  sY  ( I) 

Y3  (?)  sY  (I+l) 

Y3<  3)  = Y( !♦?) 

SAl/OR*AhS3 
SAtfPJsA  fS 

3049  CONTINUE 
3051  CONTINUE 

PRINT  93,  MSA V 
PRINT  9€,X3( 1) , Y 3( 1) 

PRINT  96 , X3(  2)  , Y 3(  2) 

PRINT  96,X3(3), Y 3 C 3 > 

PRINT  94,  SAVPJ,  SAYOR,  SAVOEL 
3061  CONTINUE 
30?1  CONTINUE 
4011  CONTINUE 

94  FORMAT  ( 1H0  ,4*X,£22.15  ) 

% FORMAT  ( 1MQ  , 6£2?«  15  ) 

9011  CONTINUE 
STOP 
FNO 
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APPENDIX  D. 


FORTRAN  LISTING  OF  THE  PROGRAM 

This  appendix  contains  the  basic  subroutine  of  this  report  which  calculates  P(A)  or  P(H)  to  3, 
6,  or  9 correct  decimal  digits. 

CALL  VALR16  (x,  y,  N,  ans,  IOP) 

where: 


j? 


I 

If 

& 


r 


W 

ft":, 


*v. 


x,  y are  input  arrays  of  the  coordinates  of  the  vertices. 
N is  the  number  of  sides  of  the  polygon.* 


f Vertic  s must  be  listed  in 


r.  See  pp.  9, 10. 


ans  identifies  the  location  where  the  Pk  is  returned. 

IOP  = 1, 2 or  3 for  3, 6 or  9 decimal  digits  of  accuracy,  respectively. 


*N  = 1 for  an  angular  region  A specified  by  three  points  given  in  counterclockwise  order  with  the 
first  point  at  the  vertex  of  A,  (See  page  12).  Note  Q < AO  v'i  2rr  for  N = 1,  but  0 < AO  < tr 
forN  > 3. 


0»-tG'Ul*-UfV»-  U M t- 


SUBROUTINE  VALR164  X,Y,N, ANS.IQP  ) 
DIMENSION  RSQ4  4) 

DIMENSION  X41),Y41),U42)»042),G(2),h42) 
DIMENSION  E45J ,E2( 10) ,E 3(15) 

DIMENSION  APH143),APH24J),APH343) 

REAL  L 

DATA  PI/3.1415  92653  56979  / 

OATA  T NOP  1/6*2831  85307  17958  / 


DATA  A LNPX/1,  1447  29885  84940  / 
OATA  Cl/.  28209  47917  73877  / 
OATA  C2/.  15915  49430  91895  / 
OATA  < Ed), 1*1,  5)  / 

,8857775185728  95E ♦SO  , 
.7593055020824856*00  , 
•6952320924352076-01  / 
DATA  (E2 (I) , 1*1,  10)  / 

,8862264700166326*00  ♦ 
.885348820003892E+00  , 

. 42 182 1 197160 099E* 00  , 
,90 5057384150 449E-01  , 

• 430 895 16098 41 38E-Q2  , 
DATA  4E341) ,1*1,  15)  / 

.8862269249314656*  08  , 

, 8 86223 7 33 16 6722 E *00  , 

• 442 85 18 9932 05 69E* 00  , 

. 1450600434030 14E* 40  , 

•3 091 9929 5521 2 10E- 01  , 
,324944543171 1 85E-Q2  , 
.1057875744806336-03  , 

.* 083355 17232 165E-06  / 
DATA  ( APH1  ( X ) , I*  1,  3 ) / 

1 2,026-7, 2,086-13, 2, 726-19  / 

OATA  < APH2CXJ,  1*1,3  ) / 

1 1,226-2,1,236-4,1.356-6  / 

DATA  4 AP  H3  ( X ) , X*  1,  3 ) / 

1 2,256-4, 2,286 -7»2,6lc-10  / 

OATA  4 RSQ  (I)  ,1=1,3  ) / 

1 6.0516,12*60605  ,19.201924  / 

NHt*H-  1 
K»1 

PHIN*0 . 

PNIKsO . 

ANS*0. 

U 41 )~X  4 2 1-X4K) 

041)*Y4  2 J-Y4K) 

IP  4 N.NE.l  I 60  TO  3131 
U42  )*X  43)  -X41 ) 

V 42 J=Y  43) -Y  4 1) 


-•9611519527780596*00 
-«  3536449806869776*00 


-, 9999507 14561036E*00 
-. 6 6 16 1 12390 433 5 7E* SO 
-.2228980556672086*00 
-,2549061 1188 A287E-01 
-.3233772396932476-03 

-.9999998997762526*00 
-•6666266705109076*01 
-.2656302  063660256*00 
-.71490 9837799889E-01 
-.1123235321484416-01 
-, 7 0426 02 433 0 90 966-03 
-.971664  86416  04616-05 
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SN*V42l*U<l)«U(2>*Yti) 
if  ( sn.&e.i*  ) i»o  ro 

ANS=-1 « 

T i3U(  1 1 
UU)=U(2) 

U (2)-T  1 

n*v«2 » 

W(2i  = V(l) 

Yll)*Ti 
SO  TO  1141 
3131  CONTINUE 

U (2) *X <N)  -XU) 

V<2)=Y  UN)-Y(i) 

3141  CONTINUE 

BGOisSQRT  i 2«*IU<l)»U(i)4W(l)*V(l>n 
BG02=SQRT(  2.*  (U(2)  *U  (2 i C 2#  (2 1 1 ) 

3 151  CONTINUE 
l*Q* 

ALAM-Q  • 

B*.5MX<K**X4KW(KI*Y(KH 
IF  ( B.GT.APHIUOP*  I GO  TO  3171 
CAPG=0. 

3161  CONTINUE 

ri*A8S<Vt2MUUI-U  (ZMfUII 
I 2*U<2l«UtU»V(2>»VU) 

PHIK=  AT  AN  2(Tt  »T2) 

ANSl=PMXK/TWOPX-CAPG 
GO  TO  3621 
3171  CONTINUE 

Gtll*U41)  ♦X1K3 ♦Will  *¥€K> 
G<2)*U(2)*X<K>+VU2»*Y(K> 
N<lJs-Y(K}*U<li*X(K>»m> 
H42**-Y(K)*U<2I*X(K)*Y(2) 

G(i)*U  UI/8G01 
G(2)»G(2)/BGD2 
H(ll»H(U/BGOl 
N12I«H|2I/BG02 

SN*U2«*<VC2»*Ulll-U(2**y<in»/4BG01*e&02i 
I IF  < SNkGT.O.  I GO  TO  3165 

CN“GU1J*G<2J*HI11*HI2) 

IF  4 C h«Gc*(U  ) GO  TC  3163 
PNlKsPI 

IF  ( G UU.IT.O*  ) GO  TO  3161 
ANSI-* 5#£RFG  1 0»H<2)  j 
GO  TO  3621 
3111  CONTINUE 

ANSI3*  5*ERFC4fl*“H(l)  ) 

GO  TO  3621 
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t^<r«v*-u>KVV  _«*«***(■ 


3163  CONTINUE 
PHIK=0  . 

ANSi=fl  . 

GO  TO  3621 

3165  IF  ( B.U.APH21I0P)  ) 50  TO  3301 

IF  { Gdi.LT.O,  ) GO  TO  3261 
IF  ( 0(2) .GE.U.  ) GO  TO  3071 
G(2)=-G (2) 

H < 2 ) = - M2) 

ALAH=P3 

IF  ( A£SIH(2))»lE«APH3(I0P)  ) GO  TO  3261 
L=.5*ERFC(0,-H(2)) 

GO  TO  3471 
3251  CONTINUE 
L = • 5 

GO  TO  3471 
3261  CONTINUE 
&U)»-G(1) 

H(l)=«H(i) 

IF  ( G (2) • LT*  0 « ‘ GO  TO  3271 
ALAHsPI 

IF  ( A as ( H ( 1 > > .Lfc.APNJ(lOP)  ) 30  TO  3251 

Us,5*£RFC(3»MU)l 

GO  TO  3**7 1 

4271  continue 

512)3-6(2) 

h(2)»-H(2) 

IF  ( ABS(H(1)) .Lt* APH31I0P)  ) GO  TO  3291 
IF  ( AdS(H(2) ) «t£t  APH3(I0P)  ) GO  TO  3281 
U = *5* (ERFC (0*  H( 1) ) -ERFG ( fl»M (2)) ) 

GO  TO  3471 
3281  CONTINUE 

U = *5M£RFCUI»M(1))-U) 

GO  TO  3471 
3291  CONTINUE 

IF  ( A8S(H(2))«LE«  APH J(IGP)  ) GO  TO  3471 
l=*5*(l.-ERFC(0,H<2))) 

50  TO  3471 
3331  CONTINUE 

C AP6-C1* (M (2)  -M(i)l“G2Mu(2)*NI2l-6(l)#Ml)  ) 
GO  TO  3161 
3471  CONTINUE 

IF  ( G*LT .kSQ(IOP)  ) GO  TO  3479 
PHINs-8. 

GO  TO  3495 
3479  CONTINUE 

IF  ( K.NE.N  ) GO  TO  3480 

IF  { FUN  *L£* lit  ) GJ  TO  3480 
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A J1*PHIN-ALAM 
60  TO  3481 
3430  CONTINUE 

SN«G<1)*H(2>-G<2)*H<1) 
CN=G(U*G<2)  + H(1)*H<2> 

A JQ-AT AN2 (SN* CN) 

PKIteAJQ 

if  ( AJO.LT.G*  ) PHIK*PI*AJQ 
3R81  CONTIN  1£ 

CAPG^Aja 
CAPH*#  5*  A JO 
*=1 
F*0. 

Ajl-M2l-HCl) 

CIRCM* AJi 

If  ( IOPtEQ.3  ) SO  TO  3631 
IF  ( I0P.EQ.2  ) SO  TO  3701 
SUM-E(N)*AJ1 
3462  CONTINUE 
H=N*i 

H(2)«H(2)«G(2) 

HU)bHU)*G(U 
T*N(2)  *H(  U 
F=F+8 

CAP0a(F*CAP6^Ti/K 
SUH*SUM*E(M)*CAPV 
IF  € N *Gc « 5 ) GO  TO  3491 

GAPG=C  JRCN 
CIRCH» CAPtf 
GO  TO  3482 
3491  CONTINUE 

ANS1»1+EXP(«I84AUNPXI I *(GAPH»SUN) 
30  TO  3621 
3495  CONTIN  t£ 

ANSl*l 

3621  CONTINUE 

IF  ( (K-NNli  ) 3631,3661,3623 

3623  CONTIN  IE 

ANSaAOSTANStABSlANSlI  ) 

RETURN 

3631  CONTINUE 
K*K+1 
KPUK*  1 

IF  ( K.NE.2  ) GO  10  3651 
ANS* AOS! ANSI) 

U<2»»X(KPD-X(K) 

V(2)»Y(KPl)*Y<K) 

PHlN*PHlN-PHIK 


5G0Z*SQRT  i 2.MU<2)»U(2)*tf(2)*tf(2))) 
80  TO  3151 
3651  CQNTINIE 
U<1J-U  (2) 

im)=v<2> 

UI2l*X4KPi)-XIK) 

C2)*YIKP1W{K> 

BG01=BG02 

BG02«SQRT(  2#*(U<2)*U<2)*¥l2)*¥t21)  ) 
GO  TO  3671 
3661  CONTINUE 
KSN 

u(i»=x*N)  *x<n 
V4D-Y  IN)  •Till 

8 GO  1=S  CRT  ( 2*  *IU(il  *U(1I  *(||  1)  (1)  I ) 

3671  CONTINUE 

PHJN*PHIN+#HIK 
ANS-ANS-AQSt ANSI) 

GO  TO  3151 
3681  CONTINUE 

SUN*E3  CM)  *6J1 
3691  CONTINUE 
H*N+1 

M42I*H  (2)  *612) 

r»H(2)«mi) 

F«F*fl 

CAP\I*CF*G  APG*TI/H 
S UN®  SUH#€  3 1 H ) *C  AP  V 
IF  I H.Gfc.15  ) GO  TO  3691 
CAPG*CIRCH 
C IRON*  CAP  V 
GO  TO  3691 
3781.  SUN*E2  (N)*AJ1 
3711  N»H*1 

M2t*HC2) *6121 
Hlt)*H(t)*GU) 

T»HI2)«HU) 

F*F*6 

CAP V*tF*CAP6*T> /N 

S UN* SUH+E  21 K) *CAPtf 

IF  I M.GE.l*  I GO  TO  3691 

CAPG-CIRCH 

CiRCN*CAF V 

GO  TO  3711 

ENO 
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